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Abstract 

It is claimed that current practices in grid convergence studies, particularly in the field 
of external aerodynamics, are flawed. The necessary conditions to properly establish grid 
convergence are presented. A theoretical model and a numerical example are used to 
demonstrate these ideas. 
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1. Introduction 

Grid convergence is an important process in verifying that discrete numerical solutions 
are valid representations of the governing partial differential equations describing the 
phenomenon under investigation. With increased emphasis on the subject of Verification 
and Validation [1], [2] and with many fluid and numerical methods journals now 
adopting policy statements on numerical accuracy, grid convergence is a topic frequently 
visited, but rarely perfonned correctly. In this note, the state-of-the-practice of 
computational fluid dynamics (CFD), particularly as it is applied for external 
aerodynamics, is reviewed and some observations are made regarding the apparent lack 
of convergence. 
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The results presented are for regular-structured grids. Similar results can be obtained 
for regular-unstructured grids (i.e. grids consisting of isosceles triangles in 2D). The 
general unstructured grid problem has not been investigated. 

2. Current Practice 

To illustrate the state-of the-practice in grid convergence, consider the drag prediction 
workshops sponsored by the Applied Aerodynamics Technical Committee of the 
American Institute of Aeronautics and Astronautics. The first of these workshops, the 
Drag Prediction Workshop-I (DPW-I), was held in June of 2001. The workshop required 
the evaluation of the force and moment coefficients by CFD of a relatively simple wing- 
body configuration known as the DLR-F4. This particular configuration was chosen 
because it had been tested in 3 wind tunnels. A total of 35 solutions based on Reynolds 
averaged Navier-Stokes (RANS) flow solvers, representing a very wide cross section of 
the CFD community, were presented at the workshop. Of these, 21 solutions were 
obtained on grids provided by the workshop organizers and 14 on grids created by the 
participants. A summary of the data presented at the workshop is given in [3] and a 
statistical analysis of the data is presented in [4]. Surprisingly, a demonstration of grid 
convergence was not a requirement of the workshop. However, in a later study Lee- 
Rausch et al. [5] assessed the grid convergence of 4 of the codes used in the workshop 
using the same DLR-F4 configuration. The second Drag Prediction Workshop-II (DPW- 
II) was held in June of 2003. The DPW-II Organizing Committee, at the request of the 
workshop participants, focused the second workshop on a more complex model [5]. The 
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new model, the DLR-F6, consisted of wing-body (WB) and wing-body nacelle-pylon 
(WBNP) configurations. The DLR-F6 wing-body configuration was very similar to the 
DLR-F4 configuration with minor changes in the airfoil sections introduced to reduce the 
chances of boundary layer separation. Wind tunnel tests for this configuration were 
conducted in 1990 at the S2MA pressurized tunnel of the ONERA center at Modane. As 
in the first workshop, standard grids were provided, but for this workshop a grid 
refinement study was required. A total of 29 data sets were submitted, however a number 
of these represented variations of the same code. The solutions presented at this 
workshop are summarized in [6] and were statistically analyzed in [7]. Additional grid 
convergence studies of the DLR-F6 configurations are presented in [8] and [9]. 

A review of the references cited above, and of the many other papers associated with 
the drag prediction workshops, reveals some curious facts. The current practice is to 
conduct convergence studies on the basis of surface integral coefficients such as the 
coefficients of drag, lift and moment. These are plotted with respect to some measure of 
the grid size, which we shall call h, raised to a power p, where p is the formal grid 
accuracy of the numerical scheme used. Finally, we notice that it is common practice in 
the literature to present grid convergence studies that do not show convergence. Hence, 
we have to conclude that the actual act of doing a grid convergence study far outweighs 
the results of the study. Figure 1 is representative of a typical figure from these studies. 
Here the moment coefficients from four of the codes used in the DPW-II are plotted 
against the grid spacing squared, since the methods used are formally second order 
accurate. 
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The practice of evaluating grid convergence with respect to h p follows from usual 
Taylor expansion studies of a numerical scheme which show that a computed quantity 
f is related to the exact value / of the solution of the differential equation by 


f e - f c + ch P + H.O.T . , 


(1-1) 


where c is a coefficient independent of h. For h sufficiently small, the higher order terms 
are negligible compared to the leading term h p and the exact solution is related to the 
computed solution simply by 


When this happens, we say that the solutions are in the asymptotic range of convergence. 
If we accept the validity of this expression and we have solutions on three grids which 
have converged in time (or iteration parameter) it is unnecessary to assume that p 
corresponds to the fonnal order of accuracy, since given three solutions Eq. (1.2) can be 
solved for the three unknowns: / , c and p. If we evaluate these quantities from the data 

reported in the DPW-II (WB case), see Table III of [10], we find that p is far from its 
formal value of two. For example, if we consider the data in Fig. (1), for the third code 
we find p= 8.54 and for the other three codes a real solution doesn’t exist, since the 
behavior of the moment coefficient is not monotonic with h. It is tempting to dismiss Eq. 
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(1.2) by suggesting that the grids used in these studies were not sufficiently fine and 
therefore the solutions are not in the asymptotic range. However, we believe that 
something more subtle is at work. 

3. The h Condition 

For two and higher dimensional problems, Eq. (1.1) is a simplification of the actual 
Taylor expansion. For a two-dimensional problem the expansion reads 


/ =/ + ah P +bh P + H.O.T. , (1.3) 

where h and h are the grid spacing in the x and y directions, respectively, and a, b are 

x y 

coefficients independent of h and h . Note that if the grid is not uniform in the physical 

space, we assume an invertible one to one mapping to a computational space that renders 
the grid spacing uniform. If that is the case, Eq. ( 1 .3) represents the expansion in the 
computational space. For both h and h sufficiently small, we have 

x y 


f ~ f + ah P + bh P . 


(1.4) 


The current practice is to replace (1.4) with (1.2) by defining 1 


1 For a three-dimensional problem h is defined as h — (h h h )* 
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(1.5) 


h = JhV . 

\ x y 


Equation (1.5) defines h to be the side of the square with area equal to the rectangle with 
sides h x , h . The definition of h is not unique. Among others, h could be defined as the 

I 2 2 

diagonal: h = Mf + h . 


There is an obvious advantage in using Eq. (1.2) instead of Eq. (1.4), since there are 
four unknowns in (1.4) and only three in (1.2), but this advantage comes with a condition. 
To be able to replace (1.4) by (1.2) it is necessary that 

h * 

-r 1 - - X for Vk, k = \, 2,3>,...,K (1.6) 

Km 

2 

where %, the grid aspect ratio, is constant over all k-grids". 

An example will illustrate this condition. Assume that we have a problem for which we 
know the exact solution and all other relevant information. Say that we know that a= 1, 
b=5,p=2, f = 1 and h is defined by ( 1 .5). If we solve the problem on three grids that 

satisfy (1.6), we find the results listed on Table I. On the table, s is the error defined by: 
s = f - f , and N x and N y are the number of grid points in the x and v directions, 

h k h k 

1 For a three-dimensional problem — — = = % forVk, k = 1,2,3, 
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respectively. If we now use the computed / k and the values of h k in Eq. ( 1 .2), we find 

p= 2, / = 1 and c=12.9, which is consistent with the original data. Now we repeat the 
exercise using three grids that do not have constant aspect ratio. The results are listed in 
Table II. If we use these results ( f c k , A ) in Eq. ( 1 .2), we find p- 2.36, f e = .999 and 

c=34.75. The results are no longer consistent with the original data. It is important to note 
that the data in Table II represents data of a second order accurate method, but because 
the grid aspect ratio is not constant Eq. (1.2) is unable to recover this fact. Note also that 
the grid points in Table II were chosen such that the ratio of successive h k ’s is constant, 

but this, despite claims in [2], does not make Eq. (1.2) valid. 

4. A Numerical Example: Ringleb Flow 

Ringleb flow describes an inviscid, compressible, transonic flow, with isentropic 
exponent of 1.4, making a 180° turn. The exact solution for this flow was found by 
Ringleb [11] in 1941. We solve the Ringleb flow problem using MacCormack’s [12] 
second order accurate, two step method on grids consisting of streamlines and iso- 
potential lines, as depicted in Figure 2. The lower boundary for the case we will be 
computing corresponds to a streamline that reaches a Mach number of 1.8 atx=0 and 
extends to the point where its Mach number is 0.6. The upper boundary corresponds to a 
streamline that reaches a Mach number of 1 .2 at x=0 and is terminated by the iso- 
potential line emanating from the lower boundary. Since the flow is symmetric about 
x=0, we only compute in the half-plane x< 0. We consider 5 grid-sets, identified by an 
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index i, for each grid-set the grid aspect ratio is constant. Each grid-set consists of 4 grids 
identified by an index k. For each grid we compute a pressure-residual to establish 
temporal convergence and an error based on the magnitude of the flow velocity defined 
by 


M N 

/(N x M), (1.7) 

m = 1 n = 1 



where q is the exact magnitude of the flow velocity, q is the computed magnitude 

of the flow velocity at grid point n, m and N, M are the total number of grid points in the 
streamline and iso-potential line directions, respectively. The results are tabulated in 
Table III. 

In Figure 3 we plot the log of svs. the log of h. The slope of the lines ^constant 
corresponds to the convergence rate p. We notice that for grid-sets 1 and 2, corresponding 
to %(.5, 1), the four grid values line-up in a straight line. Grid-set 3, %=2, deviates slightly 
from a straight line, while the deviation is more pronounced for grid-sets 4 and 5, %(3, 4). 
The significance of this is that as the grid aspect ratio increases a finer grid is needed to 
enter the asymptotic range. Thus, grid /= 1 , k= 1 is in the asymptotic range, but grid i= 4, 
A=1 is not. The last column of Table III tabulates the convergence rate p for each grid-set 
based on the slope of the two finest grids on each grid-set. The value of p rounds-up to 
2.2 for the first 3 grid-sets then slightly increases for the last two grid-sets, an indication 
that finer grids are needed to enter the asymptotic range as % increases. This is shown on 
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Figure 4 where the convergence rate is plotted as a function of h. For this figure p is 
obtained from 


P = (log(^) - log( ff i _ 1 ))/(log( A , : ) - lo g(Vi))» C 1 - 8 ) 

and h follows from 

H\ + *J /2 - ( 19 > 

Calculations with finer grids than those listed on Table III were used to create this figure. 
We also notice in Figure 3 a folding-over of grid-sets 3 and 4. This is an indication that 
with //-fixed there is a % that minimizes the error. This is depicted on Figure 5 where % is 
unfolded by plotting the error vs. the grid aspect ratio while keeping h constant at 0.0408. 
For this problem, the optimum % has a value approximately equal to 2.25. Figure 6 shows 
the error £as a function of //with p= 2.2. Only the three finest grid (k= 2,3,4) values are 
used in this figure and the error is extrapolated for /? — > 0 . These results show that Eq. 
(1.2) is valid when condition (1.6) is satisfied. What happens if this is not the case? 
Consider the following grid-set (/,&)=( 1,1), (3,3) and (4,3). This set consists of grids with 
different aspect ratios, but the error is monotonic with h. Although the data from this set 
corresponds to solutions with convergence rate of approximately 2.2, when we compute 
the convergence rate solving Eq. (1.2) we find p= 9.94. If the grids do not satisfy 
condition) 1.6), then Eq. (1.2) cannot be used to establish convergence. However, while it 
is convenient to use Eq. (1.2) it is not necessary to do so. Using the values from grids 
(i,k)=(\ ,4), (2,4) and (3,4), which clearly do not satisfy (1.6), we find from Eq. (1.4) that 
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p=2A. For this grid-set the error is not mono tonic with h and if p is evaluated by solving 
Eq. ( 1 .2), we find that p is not real. 

5. Concluding Remarks 

Today grid convergence studies with grids in the tens of million grid points are carried 
out to assess and verify the validity of numerical solutions for flows over aircraft 
configurations. In the very near future grids in the hundreds of million grid points will be 
used. These studies are expensive, time consuming and impose large demands on our 
computer infrastructure. Many of these previous grid convergence studies have been 
flawed, most probably not because of problems with the codes, but because the 
procedures used to establish convergence have been flawed. A grid convergence study 
has two primary purposes: 1) to establish convergence, and 2) to determine the rate of 
convergence. The first is significantly more important than the second. As part of a 
convergence study, we need to detennine that the solutions are in the asymptotic range 
and we need to find what the actual convergence rate is. Neither should be assumed as 
part of the process. In establishing the actual convergence rate, it is important to 
remember that the formal rate is established on the basis of a numerical algorithm that is 
used for the calculation of the field points. The algorithm is usually modified at the 
boundaries and boundary data alone might not reflect the formal rate. Finally, if we want 
to take advantage of the simplifications inherent in ( 1 .2), it is necessary to insure that the 
grids used conform to condition (1.6). 


3 This was not the case in our example with Ringleb flow. Using only boundary data we obtained 
convergence rates close to the rates we found using the field data. However, much care was taken in the 
implementation of the boundary conditions which resulted in a formally second order method at the 
boundaries. 
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Figure 2. Typical grid used in numerical calculations of Ringleb flow. 
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k 

N x 


X 

h 

e 

f c 

i 

25 

10 

2.5 

.06324 

.0516 

.9484 

2 

50 

20 

2.5 

.03162 

.0129 

.9871 

3 

75 

30 

2.5 

.02108 

.0057 

.9943 


Table I. Results from a simulation with constant grid aspect ratio. 
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k 

N x 


X 

h 

£ 

fc 

i 

25 

10 

2.5 

.06324 

.0516 

.9484 

2 

40 

20 

2.0 

.03535 

.0131 

.9869 

3 

64 

40 

1.6 

.01976 

.0034 

.9966 


Table II. Results from a simulation with varying grid aspect ratio. 
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f 

k 

N 

M 

X = e 

N/M 

pressure 

residual 

P 

1 

i 

15 

30 

.5 

0.0029603 

6.55e-007 

2.1815 


2 

20 

40 

.5 

0.0015623 

9.44e-008 



3 

25 

50 

.5 

0.0009561 

3.06e-008 



4 

30 

60 

.5 

0.0006423 

1.47e-008 


2 

1 

16 

16 

1 

0.0024823 

1.65e-010 

2.1905 


2 

21 

21 

1 

0.0013311 

1.91e-007 



3 

31 

31 

1 

0.0005612 

2.14e-008 



4 

41 

41 

1 

0.0003042 

3.80e-009 


3 

1 

20 

10 

2 

0.0017611 

1.98e-008 

2.2330 


2 

30 

15 

2 

0.0006655 

3.48e-008 



3 

40 

20 

2 

0.0003449 

1.68e-008 



4 

50 

25 

2 

0.0002095 

3.94e-009 


4 

1 

30 

10 

3 

0.0012655 

1.46e-008 

2.3558 


2 

45 

15 

3 

0.0004158 

5.68e-010 



3 

60 

20 

3 

0.0002072 

3.72e-011 



4 

75 

25 

3 

0.0001225 

2.69e-011 


5 

1 

40 

10 

4 

0.0012033 

1.45e-010 

2.4573 


2 

60 

15 

4 

0.0003636 

5.55e-011 



3 

80 

20 

4 

0.0001768 

1.02e-011 



4 

100 

25 

4 

0.0001022 

1.00e-011 



Table III. Cases computed for Ringleb flow. 
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